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Abstract 

A well-balanced scheme for a gravitational hydrodynamic system is defined as a scheme which could precisely 
preserve a hydrostatic isothermal solution. In this paper, we will construct a well-balanced gas-kinetic symplecticity- 
preserving BGK (SP-BGK) scheme. In order to develop such a scheme, we model the gravitational potential as a 
piecewise step function with a potential jump at the cell interface. At the same time, the Liouville's theorem and 
symplecticity preserving property of a Hamiltonian flow have been used in the description of particles penetration, 
reflection, and deformation through a potential barrier. The use of the symplecticity preserving property for a Hamil- 
tonian flow is crucial in the evaluation of the high-order moments of a gas distribution function when crossing through 
a potential jump. As far as we know, the SP-BGK method is the first shock capturing Navier-Stokes flow solver with 
well-balanced property for a gravitational hydrodynamic system. A few theorems will be proved for this scheme, 
which include the necessity to use an exact Maxwellian for keeping the hydrostatic state, the total mass and energy 
(the sum of kinetic, thermal, and gravitational ones) conservation, and the well-balanced property to keep a hydro- 
static state during particle transport and collision processes. Many numerical examples will be presented to validate 
the SP-BGK scheme. 

Key Words: gas-kinetic scheme, hydrodynamic equations, gravitational potential, symplecticity preserving, well- 
balanced scheme. 



1. Introduction 

Generally, flow equations with source terms can be written as 

U, + V ■F{U)^S, (1) 

where U is the vector of conservative flow variables with corresponding fluxes F{U) and S is the source term. For a 
gas flow under an external time-independent gravitational field, there exists a special solution, i.e., the hydrostatic or 
well-balanced equilibrium solution with a constant temperature and zero fluid velocity. This solution is an intrinsic 
solution due to the balance between the flux gradient and source term, i.e., 

V ■ F{U) = S. (2) 

In order to capture the physical solution for a slowly evolving gravitational hydrodynamic system, the numerical 
scheme has to be a well-balanced one in keeping the hydrostatic solution in the special situation, and has the shock 
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capturing property in the general case. Theoretically, it seems that to design a well-balanced shock capturing scheme 
for the gravitational hydrodynamic system is much more difficult than that for the shallow water equations. 

There have been many attempts to construct well-balanced gas dynamic codes which preserve the hydrostatic 
solution (lH , 15 , The schemes in ii , f5 , 2l are designed based on the condition Eq.(|2|l, such as to explicitly enforce 
this balance even for the updated non-hydrostatic solution, then use the re-balanced quantities in the evaluation of 
fluxes in the next time step. However, for a transient flow, the use of Eq.© directly in the design of the numerical 
scheme may be problematic, because in general case Eq.© is not satisfied in a physical evolution process, especially 
for flow around discontinuities. So, our aim of this paper is to design a scheme with correct particle transport and 
collision across a potential barrier, which will automatically becomes a well-balanced one when the solution is settling 
down to the hydrostatic one. But, the scheme is still accurate in capturing any general gas evolution process. 

In the past years, a gas-kinetic BGK scheme has been successfully developed for compressible Euler and Navier- 
Stokes equations without gravitational field ( ll llll2ll ). The main part of the BGK scheme is to find a gas distribution 
function / at a cell interface. Physically, the inclusion of gravitational efifect is only to change the particle trajectory. 
Therefore, it should have no much difficulty for the gas-kinetic scheme to include the gravitational effect in the modifi- 
cation of the time evolution of a gas distribution function through the particle acceleration and deceleration processes. 
Along this fine, the gas kinetic scheme (GKS) has been extended to a gravitational system ifToll . which much improved 
the solution in comparison with operator splitting method. However, mathematically, the use of a piecewise linear 
gravitational potential makes the exact solution complicated and a simplification of the numerical scheme in 1.10.1 can 
not keep a precise well-balanced solution. Therefore, the scheme presented in ifl^ is not a well-balanced one. 

In this paper, in order to design a precise well-balanced scheme we are going to approximate the gravitational 
potential as a piecewise constant function inside each cell with a potential jump at the cell interface. The detailed 
particle transport process across a potential barrier will be followed. In the construction of such a scheme, the use 
of the symplecticity property of a Hamiltonian flow and the Liouville's theorem becomes important in the correct 
description of particle penetration, reflection, and deformation processes across a potential barrier. In a previous paper 
ifTil] . following the approach of Perthame and Simeoni for the shallow water equations a well-balanced kinetic 
flux vector splitting scheme for gravitational Euler equations has been developed. However, in the above approach, 
only a few simple moments of a gas distribution function are needed, and these simple moments can be intuitively 
guessed instead of derived with a solid physical and mathematical foundation. In order to extend the above scheme to 
high-order accuracy and to solve the gravitational NS equations, a gas-kinetic BGK model with both particle transport 
and collision has to be solved. In designing such a scheme, much more high-order moments of a gas distribution 
function have to be evaluated after the interaction with a potential barrier. It becomes much harder to construct them 
intuitively. Furthermore, to model the particle transport plus collision processes through a potential barrier is much 
more challenging than that in the collision-less case. For example, around a potential jump at a cell interface, a multiple 
equilibrium states have to be constructed on both sides of a jump. In the construction of such an equilibrium state for 
the BGK model, the second law of thermodynamics has to be satisfied. 

The paper is organized as follows. In section |2] we will present the basic physical principles about the particle 
interaction with a potential barrier. The symplectic principle plays an important role in the design of the well-balanced 
scheme. Section[3]gives a brief review of the previous BGK scheme without external forcing field. Section |4]presents 
particle transport mechanism and the construction of a symplecticity preserving BGK for the gravitational gas dynamic 
system. Section|5]is about the theoretical analysis of the schemes, such as the necessity of using an exact Maxwellian 
and the well-balanced property. Section|6]shows the numerical tests. The last section is the conclusion. 



2. Particle transport mechanism across a potential barrier 

In this paper, the gravitational potential (p is modeled as a piecewise constant function. With 0y in yf/i-cell and 
(pj+\ in ij + \)th cell, there exists a potential jump at the cell interface, i.e., A0y+i/2 = (pj+\ - (pj. Now what we need 
to figure out is the eff'ect on an initial gas distribution function next to the potential barrier when the particles move 
towards the barrier The associated physical process could be reflection or penetration of the particles from the barrier. 
What we have to evaluate is the relationship between the moments of the gas distribution functions before and after 
interaction with the potential barrier. Since all particles are located next to the potential jump, the modification of 
the particle distribution function happens instantly. Therefore, once a time-dependent gas distribution function next 
to the potential barrier is given, the corresponding distribution after particle collision with the potential barrier can be 
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evaluated at that moment. Since the potential jump only affects normal velocity and its moments, so in this section we 
only consider distribution functions with 1-D velocity. The results obtained in this section will be used in this paper 
many times on the construction of symplecticity-preserving scheme. 

For an initial gas distribution function f(u) next to a potential barrier and these particles impacted with the potential 
jump, the particle velocity u changes to u', and the distribution function becomes /(m')- We are going to use the 
following three physical principles to find the relation between the velocity moments of /(m') and /(m). 

a. Hamiltonian preserving property: the Hamiltonian function Hofa particle keeps a constant, where 

H^^u^ + <p{x). (3) 

This is actually the energy conservation for a particle movement under a conservative potential field. Since we only 
consider the interaction of a particle with a potential barrier at an instant of time, there are no collisions between 
particles. Therefore, the energy conservation for individual particle is precisely conserved, i.e., 

i«2 + ^ ^ i(„')2 + (4) 
from which the relation between u and u' can be obtained. 

b. Liouville's theorem: the probability density of a particle in phase space keeps a constant along its movement 
trajectory, 

/("') = /(«). (5) 

In other words, the particle isn't lost or created during its impact with the potential. 

c. The symplecticity preserving property: for a Hamiltonian phase flow, we have 



^ J" dx'du' ~ J" J" dxdu, 



(6) 



where D' and D are the phase volume on the trajectory of the Hamiltonian phase flow. 

During the impact of the particles with the potential barrier, we can specially choose D - {ui, M2) x (uti,ut2), then 
D' = (m'p Mj) X (u'ti, u't2) since D and D' are on the trajectory of the same particle. Therefore, Eq.® goes to 



u'du' = 



«2 

udu. (7) 



This relationship will be the most important one in the construction of the moments between between f{u') and 
/(m). Therefore, the developed scheme in the present paper which uses this relationship will be called symplecticity- 
preserving scheme. 

With the above three physical principles, we can derive the relationship between the nf/j-order velocity moments 
of f{u') and that of f{u). From (|5]l and O, we have 



I f(u')u'du' - I f(u)udu 

J u'. Jui 



(8) 



Moreover, Q tells us that u' is a function of u, i.e., u' - u'(u). So, combining with ([8]), we can get a general 
formulation. 



nth-order M moment = I f(u'){u')"du' - I f(u)(u'(u))"^^udu. 



(9) 



which connects the moments of the distribution functions before and after impacting with a potential barrier at an 
instant of time. The above distribution function can represent the portion of particles which are reflected or penetrated 
at the barrier 
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3. A review of gas-kinetic BGK-NS scheme without external forcing field 



The BGK equation without external forcing field in 2-D is 



g-f 



(10) 



where / is the gas distribution function and g is the equilibrium state approached by /, V/ is the gradient of / with 
respect to x, x = (x, y), and u = (u, v) is the particle velocity. The particle collision time r is related to the viscosity 
and heat conduction coefficients, i.e., t - fi/p where jj is the dynamic viscosity coefficient and p is the pressure. The 
relation between mass p, momentum (pU, pV), and energy pE densities with the distribution function / is 



( P\ 
pU 

pV 

■pE. 



where 



(11) 



- d^id^2---d^K, and K is the number of degrees of internal freedom, i.e., K - {A - 2y)/('y - 1) for 2-D flow. Since 
mass, momentum, and energy are conserved during particle collisions, / and g satisfy the conservation constraint. 



///< 



(g - f)4fadudvd^ = 0, a =1,2,3,4 

at any point in space and time. The integral solution of ( fTOt is 

1 r' 

fix, t,u,^)^ - g(x', f, u, ^)e-'-'-''^'^dt' + e-"^fo(x- ut, u, ^), 

T Jo 



(12) 



(13) 



where x' - x - u(t- t') is the particle trajectory. The solution / in ( fT3] l solely depends on the modeling of /o and g. 

For a finite volume scheme, we need to evaluate the fluxes across a cell interface in order to update the cell averaged 
conservative flow variables. In the BGK scheme, the fluxes are defined by 



F \ 
Fpu 
Fpv 

\FpEJ 



///" 



4tfdudvd^, 



(14) 



which depends on the gas distribution function / in Eq.(fT3]) at the cell interface. Let's consider the construction of the 
distribution function at the cell interface -?,+i/2 = (j^j+i/a,}'/), where Xj+1/2 is the location of the cell interface center in 
the physical domain. Locally, around this cell interface, with the assumption of the x-direction as the normal direction 
and y-direction as the tangential direction, based on the BGK model a solution in this local coordinate can be obtained. 

By using the MUSCL-type limiter, a discontinuous reconstruction of the macroscopic flow variables can be ob- 
tained around the cell interface (see fig[T]|. The initial gas distribution function /o in (fT3]) on both sides of a cell 
interface can be constructed as 



fl^{x, u,^) = g(,(l + a'(x - Xj+1/2) + b'iy - y,) - T(a'u + b'v + A')), x < Xj+1/2, 
f[^{x, it, ^) = g|)(l -I- a'ix - Xj+112) + b'iy - yi) - T(a''u + b' v + A')), x > Xj+1/2, 



(15) 



where the Chapman-Enskog expansion up to the Navier-Stokes order has been used in the above initial reconstruction. 
Here g'^ and g'^ are the corresponding Maxwellians to W' = (pi, ipU)i, (pV)i, (pE)i) and W - (pr, (pU)r, ipV)r, (pE)r) 
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at both sides of the interface. The Maxwellian distribution function corresponding to W - (p, (pU), (pV), (pE)) has 
the form 



K+2 

g^p(fj^' e''««-^>'+(-^>'+f'', (16) 

where A is equal to m/2kT, m is the molecular mass, k is the Boltzmann constant, and T is the temperature. The 
equilibrium distribution functions around the cell interface can be modeled as 

I ^ ^ I —I ~' 

g'(x,t,u,^) = §^+1/2(1 + fl (x - xj+1/2) + b (y - yd + A t), x < Xj+1/2, 

(17) 

/(f,f,M,^) = ^^+1/2(1 + a' (x - Xj+1/2) + b iy-yd+A t), x > Xj+1/2. 

In the case without external forcing term, ^'+1/2 ™d g'j+1/2 ™ the above equation are the same distribution functions, 
i.e., g'j_f_i/2 - g'j+i/2 ^^^^ which can be obtained using the conservation constraint (fT2l l at ic- -?/+i/2 and f — > 0, 

(18) 

= //X>o /()<^^J+i/2' «' ^Wudvd^ + ///^^^ /o (^;+i/2, M, ^)(AdMdvd^. 

Therefore, at the cell interface the final distribution function can be fully determined using the integral solution ( fTlT l. 
The final distribution function can be written as 

fiXj+l/2,t,U,0 

( f{Xj+ll2,t,U,^) M > 0, 

\ rixj+ii2,t,u,^) u<0, (19) 

, 7 X' g'{xj^ii2 - u(t - n, f, a, ^)e-('-*'>l'dt' + e-'/Vo (-^i+i/2 - ut), u < 0, 
which can be used to evaluate the fluxes 

(20) 

= //X>o uf'(xj+ii2, t, M, ^)iffdudvd^ + uf{xj+ii2, t, u, ^^|/dudvd^. 

The update of the cell averaged conservative variables becomes 



where F'j_yi^{t) ... F''.^yi^{t) are the fluxes at the center of the cell interfaces. 

The definitions and constructions of all parameters related to the spatial and temporal slopes, such as a, b and A, 
can be found in 1 11] and 1 12]. 

In summary, at the cell interface X;+i/2 we can construct the equilibrium distribution functions g^j^ii2 ^"d g'j^ii2 
from initial distribution and /q. Also, we can find fluxes F'j^-^^^i^) from the integral solution /' and 

f. Without external forcing field, all the particles running into the cell interface can freely cross it. Therefore, 
the equilibrium states and fluxes at the interface have unique values, i.e., §'_|_i^2 ~ S'j+1/2 ^j+i/2*-^^ ~ ^j+i/2*-^^- 
However, with the approximation of constant potential inside each cell and a potential jump at the cell interface, the 
modeling of equilibrium state g around a cell interface has to be considered separately on different sides of the cell 
interface, where g'j_f_i/2 S'j-^i/2 general case. But, the mathematical formulae described in ( fTTb and the integral 
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solution in Eq.(fT9]l can be still used. One of the main reason for the validity of the integral solution is that there is no 
gravitational force inside each cell. However, the construction of the equilibrium states and the calculation of fluxes 
will not be as simple as that in ( fTSl ) and (l20l i. In the evaluation of the equiUbrium states and the fluxes, the physical 
principles for the particle transport discussed in the last section have to be used. In the next section, the determination 
of g and fluxes will be described. 

4. The symplecticity preserving BGK(SP-BGK) scheme 

In this section, we will construct a well-balanced gas-kinetic scheme for hydrodynamic equations under gravita- 
tional field. In order to clarify the concepts, we are going to use a similar procedure as that of the construction of the 
BGK-NS scheme without external forcing field. 

4.1. The initial data reconstruction 

For a hydrostatic solution, the flow variables satisfy the conditions, 

U — 0, V — 0, A — constant, Ba — constant, (22) 

where Ba - pe^'**. In order to avoid introducing errors in the initial reconstruction for the hydrostatic case, it is 
reasonable to use the variables {U, V, A, Ba) in the reconstruction. More specifically, we firstly apply a MUSCL-type 
limiter to reconstruct the slopes of {U, V, A, Ba), i.e., {S u, S v, S a, S Ba) inside each cell. Since 

Ba 1 , . K + 2 

P-^,'PE-2P^U' + V') + ^p, 

we can get the corresponding slopes for other flow variables, 

Sp = —^Sbo - '2p4>SA,Spu = SpU +pSu,Spv = SpV + pSv, 



1 



-(U^ + V^) + 



K + 2' 
~AA 



Sp +p 



USu + VS 



V ■ 



K + 2 



where {S p, Spu, Spv, SpE) are the slopes of (p, pU, pV, pE) inside that cell. Therefore, we can reconstruct (p, pU, pV, pE) 
in each cell using their cell averaged quantities and the above slopes. Here, all slopes become zeros when the initial 
flow is in a hydrostatic state, and the reconstruction will not introduce numerical errors. In the general case, the above 
reconstruction works as well. 



4.2. The gas-kinetic SP-BGK scheme 

With the modeling of piecewise constant gravitational potential inside each cell, i.e., (pj inside the jth cell, there 
is a potential jump at the cell interface ij+i/2- It is obvious that the distribution function / also satisfies the equation 
( [Tol l inside each cell since there is no external forcing term inside each cell. Therefore, the similar framework used 
in the constructing BGK-NS scheme can be extended here to design the SP-BGK scheme with gravitational field. 
For example, with the initial reconstruction, the non-equilibrium states around each cell interface can be obtained. 
Also, due to the potential jump, the equilibrium states are different in the left and right hand sides of the interface, 
but the integral solution of the BGK model can be still used in the construction of the local solution separately around 
the cell interface. However, at the cell interface, we have to consider the effect of the potential jump on the particle 
movement. Since the equilibrium states, g'j^ii2 ™d •?j+i/2' '■^^ fluxes, F'j^-^^^i^) ^"d F'j^-^^^it), involve the particle 
interaction with the potential jump, we will show that g'j+1/2 ^ ^7+1/2 ™ Eq.(fT7b(see figlH, and F'j^-^^^i^) ^''j+iji^^^ ™ 
the general case. Their determination depends on the particle transport modeling. The potential jump gives a critical 
speed Uc - '<J'2\<pj - 4>j+i\^ which provides a threshold for the particle movement. Because of the potential jump, not 
all particles running into the cell interface could go through freely. Some may be reflected due to less kinetic energy to 
overcome the potential barrier (see figO). For these particles passing through the cell interface, their momentum and 
energy need to be modified due to particle acceleration during the transport process. 
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Without losing generality, we only discuss the case of 4>j < 4>j+i in this subsection. Using similar methods and 
ideas, all the formulae for the case <pj > (pj+i can be easily obtained. Let's assume the initial reconstructed gas 
distribution at a cell interface before the interaction with the potential jump is 

f fj(xj+U2,t,u,0^ u>0, 
/(f,+i/2,f,i?,f) = j ^ ^ (23) 

Starting from the above distribution function, the particle collision with the potential jump changes distribution 
functions to "'^^ ^'^'^ "'^-^ ^^^^ ^'^'^ right hand sides of the cell interface respectively, which can 

be represented as 

' fj(Xj+l/2,t,U,0y M > 0, 

t, u, ^), 0>u> -Uc, (24) 

[ fj+l(Xj+\/2'f''^'0' u<-Uc, 

and 



[ fj(Xj+l/2,t,U,^), M > 0, 

/J,i/2(f,i?,^)= ^ ^ (25) 

I fj+\iXj+l/2J,U,^), M < 0. 

The definition of the above distribution functions is from the following physical consideration (see figO. Because 
the potential jump is only at the normal direction of the cell interface, it only affects the normal particle velocity, u. 
In (l24l) . fj is the distribution function of the reflected particle in the jth cell with the original distribution function 
fj which has a positive particle velocity less than Uc- Here fj_^i is the distribution function of the particle in the jth 
cell coming from the {j + \)th cell with the original distribution function fj+\ with negative particle velocity. This 
particle has been accelerated in the negative normal direction after passing through the cell interface. Also, f j is the 
distribution function of the particle in the {j + \)th cell coming from the jth cell with the original distribution function 
fj and positive velocity higher than Uc- This particle has been be decelerated in the positive normal direction after 
passing through the cell interface. Therefore, the effect of the potential jump modifies the distribution function, but the 
particle velocity moments of the modified distribution function and the original ones are related through the physical 
principles which have been introduced in section |2] 

Here, we will show the procedure of the SP-BGK scheme first, then clarify the detailed derivation of the formulae 
for equilibrium states and fluxes. 

Using particle free transport mechanism in Eg.dTSTl for the initial gas distribution function /o, i.e., /,(x/+i/2, f, u,^) - 
/y(xj+i/2 - Mf, M, ^) and fj+\{xj+\i2, t, u,^) = ,fl^{xj+ii2 - ut, u,^), and due to their interaction with the potential jump, 
the initial condition will be changed according to Eq. (l24l) and (1251 ). from which two sets of conservative variables at 
different sides of the cell interface can be obtained. 

Km = ///I4i/2(f = 0,<7,^)^d«dv^ 

= // fo°° M^nm^ f = 0' «' ^)'AdM + // fj(xj^i/2, t^Oj, ^Wudvd^ (26) 

and 

Km = ///I/;;i/2(' - o,«-,^)^d«dv^ 

(27) 

= JJ£°° fj(xj+\/2, f = 0, M, ^)i//dudvd^ + //£^ fj+\{xj+\i2, t = 0,u, ^)i//dudvd^, 
from which, two Maxwellians g'j_f_i/2 ^^'^ •?j+i/2 ™ equilibrium states (fTTl l can be fully determined. Then, following 



the method used in the development of BGK-NS scheme 01211 , the final gas distribution at the left and right hand sides of 
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a cell interface, i.e., /' and f in ( fT9] l. can be obtained. When choosing the integral solutions as the original distribution 
functions, i.e., fj{xj+i/2,t,u,^) - /'(x,+i/2, f, «, ^) and fj+i{xj+i/2,t,u,^) - f''(xj+i/2,t,u,^), and considering their 
interactions with the potential jump, these distribution functions will be modified as Eq. (l24l) and dZSl ). from which the 
corresponding fluxes at different sides of the cell interface can be determined, 

^j.i/2(o = //rr <i/2(^'«'^)'Ad«dvd^ 

Jf^°° ufj(xj+i/2, t, u, 4)ipdu + J j*^ ufj(xj+i/2, t, u, ^)if/dudvd^ (28) 



and 



+ IILco ' «/j+i(^>+i/2: ^ u, ^)(AdMdvd^, 
^;+i/2(0 - nC ^f;.U2(f^ «,^)^d«dvd^ 

= //X^°° "7j(-^;+i/2' f> «> ^)i/'dMdvd^ + jj£°^ ufj+i(xj+i/2, t, it, ^)iffdudvd^. 



(29) 



Note that due to the potential jump, in general we have / j ^ /2 ^'^d F\+\i2 ^'m iv Finally, we can use (|2TI) to 
update the cell averaged conservative variables. 

In the above formulae ( l26l l. (l27T i. (l28T l and (|29T l, we need to find the «f/i order velocity moments of the modified 
distribution functions, fj, /^^j and which can be evaluated from the moments of the original distribution funcions 
fj, fj+\ and fj respectively by (|9]i- Let's figure out how to evaluate the nth order normal velocity moments of /,(«), 
7/+i(m) and/^M). 

a. The nf/z-order normal velocity moments of // 

Recall that fj is the distribution function of the reflected particle in the jth cell. Assume that the normal particle 
velocity is u before the reflection, and the distribution of the particle before reflection is fjiu) with Q < u < Uc- After 
the reflection, its velocity becomes u' and u' - -u, for these particles, ^ gives 

fj{u')(u')"du' = I fj(uM-u)"-^du = I //M)(-l)"M"dM. (30) 
Uc Juc Jo 

b. The nf/i-order normal velocity moments of /^^i 

f is the distribution function of the particle in the jth cell coming from the (j+l)th cell. Its distribution function 
before crossing the potential jump is fj+i with normal velocity u < 0. After passing through the interface, the normal 
velocity changes from u to m', where u and u' are related by the Hamiltonian preserving property, i.e.. 



1 2 lw^2 

-u + (f>j+i ^ 2^u) + I 



bo, M gives 

fj^i(u')(urdu' ^ fj^du)(-ir-'u(u^ + U^J"-'^'^du. (31) 

oo %J — oo 

c. The nf/z-order normal velocity moments of fj 

f I is the distribution function of the particle in the (j + l)th cell coming from the jth cell. Its distribution function 
before passing through the potential jump is fj with normal velocity u > Uc- After passing through the cell interface, 
the normal velocity changes to u'. The relation between u and u' becomes 

1 2 1. ,^2 

-U +4>j = -(m) +(pj+l. 



So, u' = V(«?^^' Eq-® deduces 

I fj(u')(u')"du= I //m)m(m2 - f/2)("-i'/2jM. (32) 

Jo J(7£ 



Based on the above moment evaluations, we can get the formulae for W'.^j^j' ^/+i/2' ^^'^ 
( I32I 1 for the case (f)j < (f>j+\- The formulae for the case 0j > (pj+x can be found similarly. All the formulae are 
given in the appendix for both 1-D and 2-D cases. Therefore, the SP-BGK scheme is presented. 

4.3. Limiting Cases 

a. The 1st order SP-BGK scheme 

When all the slopes in the reconstruction are zeros, and all slopes a, b and A of the distribution function in ( fTSl ) and 
(fTTl l become zeros, the SP-BGK scheme becomes a 1st order scheme. Now, the distribution function in ( fTSl) becomes 



/(x,+ l/2,f,M,^) 



+ M<0. 

Or, with the definition of a small parameter e, i.e., < e < 1, the distribution function becomes 

(1 -«)^;+l/2 + '?^0' 



fiXj+l/2,t,U,0 



which is called the Ist-order SP-BGK scheme. 



(33) 

(1 -e)^;+i/2 + fi^J)' «<0' 



b. The SP-KFVS scheme 

When the collision time r goes to +00, the distribution function in ( fT9] l becomes 



f'(Xj+i/2,t, 




u>0. 


f''{Xj+l/2,t, 




u <0, 




- lit). 


M > 0, 


( fo(^j+l/2 


- ut). 


M < 0. 



(34) 



The above solution solely comes from free transport and there is no contribution of the equilibrium states g in the 
integral solution /. It equals to solve 

/, + i?-V/ = 

directly when the initial distribution function is modeled as ( fTsT l. In other words, we don't consider particle collision 
here, and needn't to model the equilibrium distribution function g in ( [TtI i. This is exactly the same scheme introduced 
in [14], which is called SP-KFVS scheme. It is actually a limiting case of the SP-BGK scheme. 

In this section, with the assumption of piecewise constant gravitational potential, a SP-BGK scheme is presented. 
As will be presented in the next section, the SP-BGK scheme is a well-balanced scheme for the gravitational hydro- 
dynamic system. This is the first well-balanced scheme, which has the shock capturing property as well in the general 
case. 



5. Theoretical analysis 

For simplicity, we are going to prove all the theorems in the 1-D case. But all the conclusions still hold for higher 
dimensions as well, because there is no dynamic difference in higher dimensions when the potential jump is modeled 
as a piecewise constant function. 
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In the current scheme, the updated flow variables inside each ceU are the mass, momentum, and energy densities 
(kinetic + thermal ones). The gravitational energy is not explicitly included. However, for an isolated gravitational 
system, the total energy (kinetic + thermal + gravitational ones) conservation is a necessary condition in order to get 
a correct physical solution. In the following theorem, we are going first to prove that the conservation of total energy 
in the current kinetic scheme is satisfied. 

Theorem 3.1: The SP-KFVS and SP-BGK schemes are mass and total energy conservative schemes. 
Proof The only difference between the SP-KFVS and SP-BGK schemes is that they have diff'erent original distri- 
bution functions fjiii) and fj+[(u). However, whatever /^(m) and fj+\{u) are, the mass and total energy are conserved 
when the fluxes are calculated by (|92] | and (l93l l or (l94l l and (l95T l in the appendix. The concept of conservation of a 
variable means that the change of that variable in any fixed domain depends only on the fluxes across the interfaces 
of that control volume. In the following proof, we assume the control volume consists of many cells between the 
cell index Ki and K2, where Ki < K2- Then, we need to prove that the change of the mass and total energy in the 
control volume depends only on the fluxes at the interfaces xki-i/i and xk2+i/2- Without losing generality, we assume 
<pj < <pj+i everywhere. 
Mass conservation: 

For mass, in each cell we have 

where F'^.'' „ are the mass fluxes. The total mass in the control volume is E'^V Pi, and 
From ( l92l i and (l93T i. we have 



(37) 

^ / //«)«dMd^ -H / j^fj+i(u)udud^ 
- f 

^7+l/2,p- 

Therefore, from (l36l l and (iJTl l. 

1.%, pT = ^Ik. p] + h V - J dr. (38) 

which gives the mass conservation in the computational domain. 
Total energy conservation: 

The kinetic energy and thermal energy, i.e., pE, is updated by 

(p£)f' = (p£)« + i- - dt, (39) 

where F'^'j^ij2pE fluxes of pE. Because the external potential is independent of time, the potential energy, i.e., 

p0 is updated by 

With the definition of total energy TE - pE + p(f>, we get 

(41) 

+^}-l/2,p£ ~ ■^'■+l/2,p£] 
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The updating of the total energy in the control volume (i.e. Tj^I^^ TEj) becomes 

y^v = y'^v te" + i- f'""' y^v If' , <*,■ 

-F' cpi + F''. ,„ rldf. 

j+\/2,p"rJ J-U2,pE j+l/2,p£j 

According to ( l92l i and (l93T l, we get 

+ f /_l/>+i(«)5(«(«' + f/') + OdMd^, 



(42) 



(43) 



+ 

A direct calculation gives 

So, from (l42l i and (l44l i. the total energy update becomes 

y'^v r£"+i = y^v te". + -rL f'""' [fl ,„ <*^, 

^j=Ki J ^J=Ki J ax Jt„ I Ki-l/2,pT'^l 

~Fk2+\I2,p^^^ Fk,-i/2,pe ~ ■^L+i/a.pfi] d^ 



(45) 



which guarantees the total energy conservation in the whole computational domain. Based on the above proof, the 
SP-BGK and SP-KFVS schemes are conservative methods. Therefore, the above two schemes can give the correct 
shock location even with the external gravitational forcing terms. This is a generalization of Lax-Wendroff theorem to 
the system with gravitational source term fl- 

Lemma 3.2: The density p{x) in a hydrostatic state under the gravitational field 0(x) satisfies 

p{x) = Cxe-^^^1'^-'\ (46) 

where C\ and A are constants. 

Proof For a hydrostatic solution under the gravitational field 0(x), we have 

Px - -p<f>x, T = constant, U -0. (47) 

Since T = constant and A - m/2kT, we know A = A, where A is also a constant. Then from ( |47] ) and the ideal gas 
equation of state 

1 

P - —P' 
2A 

we have 

1 

T^Px = -pif>x- 
LA 

Therefore, with a constant, C\, the solution becomes 

p(x) = Cie-^^-^W. 

Remark: without losing generality, in the following proofs, we let C\ - \ for the hydrostatic solution. So, in the 
hydrostatic case, the state has the form 

p = e-2-''^'W, = 0, (48) 
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where 1 is a constant. Numerically, if we let the potential <^(x) be a constant, (pj, in the jth cell, then 
where pj and Uj are cell average quantities in that cell. 

Lemma 3.3: For the two equilibrium states Wj^j^j = ^+i/2' (/'^);+i/2' (P^);+i/2) = ^'^>+i/2' (P^)>+i/2' (/'^)y+i/2)' 

they have the following properties when the initial flow is in a hydrostatic state. 

1. Both velocities are equal to zero, i.e., 

Km = Km = 0- (50) 

2. They have the same temperature at both sides of all cell interfaces, i.e. 

^j+m - "^j+m 



'^',■+1/2-4+1/2-'^' (51) 



where A satisfies 

p£-ipf/2=p^, (52) 

macroscopically with K - (3 - y)/{j - 1) in 1-D, and A has the constant value A of the hydrostatic solution. 

3. The densities at the same cell interface satisfy 

4. In the same cell, 

p'nm=Pj-m (54) 

Proof As the definition, W'.^j^j ^^'^ ^j+i/i determined by dSST l and ( [89] l or ( |90] l and ( |9TT l for < ^y+i or 
<f)j > (pj+i when fj{u) - gj{u), where gj{u) is a Maxwellian corresponding to the cell average conservative variables, 
{pj, (pU)j, (pE)j). Here, we only prove the case for (pj < ipj+i- The other case can be proved similarly. From direct 
calculation, we can get 

p ' A r*^ A A /^+"' j 

P/+i/2 = T+^/(-)- e-'"-du-pj^i(-r'U,+pj^a(-r- e-^^'Jt+U^dt, (55) 
2 ■ n J_u^ 7T Jo ' 

PJ.U2 - pA(^)^- £7 A/^df + (56) 
(pf/);.,i/2 = ipUYj.ui = 0, (57) 

and 



(58) 



where f/^ = -y/2((^;+i - 0;). 



(59) 
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1. From i55[ and ( |56] |. we can easily see that p'+i/2 ^ P'j+i/2 ^ whenpj > andp^+i > 0. Since U = pU/p, 
from (ISTT i. we know that 

f/;.I/2 = f/;.l/2=0. 

2. From (|55]l, 

-n' + + PJ*' [In + P' fl r° ^-''"'h« 



„/ K+l _ g 



Since (p£)^.^j^2 " 5P^>i/2(f^i>i/2)' = P^-.i/2 if^ t^>.i/2 = 0' have 

1 



Therefore, substitute ( l49l l, dSSt and ( l60b into dHTT l, we get 

~n I 1 J . Pj Pi [^„-~AUljr 



Because e ' is a mono tonic increasing function on [-f/c, 0], so 



Then we know the summation in the brace {...) of ( l62b is strictly larger than zero. Therefore, 



4+1/2 - ^■ 



has to be satisfied. 

Similarly, we can have 



Again, the summation in the brace {...) is strictly larger than zero. So, 



3. It is easy to prove that 



and 



So n'' — n' ^-2i(0,+i-0j) 



XO _ ^U; 
(7( Jo 



(60) 



(p4.i/2-p'-.i/2T7— =0. (61) 
7+1/2 



(62) 



^ f e-'"-du-^J^e-'"^U,>0. (63) 



^dM, (64) 



(65) 
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lefl:.x=t-U; ■,right:x=t+Ul _|_(, 



Jo ' J-(/< 

Therefore, from (|65]|. we can conclude that 



Pj+1/2 - Pj+l/2^ 

<s^> Pj J^^ e--^"'du - pje-'^"' Uc + PjA £r e^'^" ^^xdx = pjA j^'" e^'^^ yfxdx, 



iXpj e''^"~u^du + PjA e^'^^ ^|xdx = pfA e'^^ y[xdx. 



From ( 165b . we know that the last equality holds. Therefore, 

P>+l/2 = P}-l/2- 

Remark: the above lemma, especially part 2, illustrates that starting from a hydrostatic state with the same temper- 
ature, the constructed equilibrium states at both sides of a cell interface have the equal temperature as well. In order 
words, in the hydrostatic case, the particle interaction with the potential barrier and the particle collisions among them- 
selves never alter the equilibrium temperature both sides of a cell interface. This is consistent with the second law of 
thermodynamics. Otherwise, the temperature differences generated by the particle collisions could be used drive an 
engine and a pure work could have been extracted from an initially isothermal system. This violates the 2nd-law of 
thermodynamic s . 

Theorem 3.4: For a well-balanced kinetic scheme, the equilibrium distribution function must be an "Exact Maxwellian". 
Proof In order to keep the hydrostatic solution (|49] | the numerical mass flux at both sides of a cell interface must be 
zero. 

Without losing generality, we only consider the case for > ^ j. Since the gas must be isotropic, we can assume 
the equilibrium distribution function is p(x)G{u^) and define a - -^2(0^+1 - <pj), then we require 



(66) 



r°° 2 r° 2 

^'j+\/2p~ I PjG{u)udu+ I pj+\G(u )udu - 0, 
where F'j^-^^^p is the mass flux at the right side of the interface. Because of ( |49] l, we have 

- I Gix)dx + e-'^"' I G(u^)udu^Q. (67) 
Take the derivative of ( |67] | with a-, we get 



2 r" 

G(a^)-Ae-'^" G(u^)udu ^ 0. (68) 



It is obvious from (|68] | that 

G(fl2; 
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G(a^) ~ e-^"', (69) 



which means that the equihbrium distribution function is an exact Maxewelhan distribution. 



Theorem 3.5: Both the Ist-order SP-KFVS and SP-BGK schemes are well-balanced schemes. 
Proof In order to prove a scheme to be a well-balanced one, we only need to verify that the scheme can keep the 
hydrostatic solution ( |48] ) forever. Numerically, the initial condition for this case is given by ( |49] l in the jth cell. At the 
next time step, the above solution must be kept by the well-balanced numerical scheme, i.e., WJ^' = W". From ( |2TI) . 
we must have 

Fj-m = F'nm- (70) 



Therefore, to complete the proof, we have to show that mass fluxes (^^JVi/2p)' momentum fluxes (^^+i/2p{/) energy 
fluxes (F'jl^^2pE^ satisfy the condition dTOl l respectively. 

The lif-order SP-KFVS scheme: the original distribution function at the cell interface is 



( gj(u), M > 0, 

/(x,+ l/2,f,M,^) = j (71) 
[ gj+l(u), M < 0, 

where gj{u) is the Maxwellian corresponding to (pj, {pU)j, (pE)j). The proof is only a direct calculation of the fluxes 
at the interface using ( |92l ) and ( |93T l or ( |94l ) and ( |95l l in two different cases for (pj < (pj+i or (f)j > (f>j+\- Also the initial 
hydrostatic condition ( |49l ) will be used. The results are the followings. 
a. For mass flux. 



b. For momentum flux. 



c. For energy flux. 



fU/2,,-F^.U2,p-0. (72) 
F'j+\/2,pU = F'j_i/2,pu = (73) 



F'j+i/2,pE - F'j+nxpE - 0- ("74) 
Hence, the first order Isf order SP-KFVS scheme is a well-balanced one. 



The 1st order SP-BGK scheme: the original distribution function is 

''■+1/2 



f{Xj+l/2,t,U,^) 



(1 - e)gj(u) + eg'.^^Ju), M > 0, 



(75) 

(1 - e)gj+i(u) + eg^._^j^2(M), M < 0, 



where e is a constant between and 1, gjiu) is the same as in the proof for the 1st order SP-KFVS scheme, gj_^i/2 

and g'j+i/2 ^'"^ '■^^ equilibrium states corresponding to W'^j^j ^^'^ respectively. Here, Wj^^j2 ^^id ^'"^ 

macroscopic variables calculated by (l88T l and ( |89] l or ( |9oi and (|9TT l when 

//■(«) = 8jM and //+i(«) = gj+i(u)- 
So, the fluxes are the linear combination of two kinds of fluxes Fi and F2 calculated by 



/i = I and /2 = 



respectively. 

From the above proof for the 1st order SP-KFVS scheme, we know that the first kind fluxes Fi can satisfy (iTOl l 
itself. Therefore, we only need to prove that F2 can satisfy (iTOt . too. Note that in the proof for the 1st order SP-KFVS 
scheme, the hydrostatic initial condition is the key. But from the Lemma 3.3, we can see that the equilibrium states 
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also satisfy the hydrostatic initial condition. So, similarly, we get the following results for the fluxes corresponding to 
/2 from a direct calculation by using ( |92] i and ( |93T l or ( |94l ) and ( |95] ) in two different cases for < (f)j+i or (pj > (pj+i- 
a. For mass flux, 

^',■.1/2.0 = F^^no = 0. (76) 



b. For momentum flux. 



(77) 



Based on Eq.dSTTl and ( l54l l. 
c. For energy flux. 



F'j+l/2.pU - F'j-i/2,pU- (7^) 



F'j+l/2,pE - F"j+l/2,pE - 0- ("79) 

From all the above proofs, we can conclude that both the Ist-order SP-KFVS and SP-BGK schemes can keep the 
initial hydrostatic solution forever Therefore, they are well-balanced schemes. 

Remark: The 2nd order SP-KFVS and SP-BGK schemes are well-balanced schemes. 

We use (U, A, pe^'^'^) to do the reconstruction. All the three variables are constants when the solution is in a 
hydrostatic state. So, the slopes are all zeros after using the MUSCL-type limiter In other words, the 2nd-order 
schemes go back to the 1 st-order method when the solution is in hydrostatic state, which can be kept forever. Therefore, 
the 2nd-order schemes are also well-balanced schemes. 



6. Numerical examples 

In this section, we will present numerical results of four 1-D examples by using 1 " and 2'"' order SP-KFVS and SP- 
BGK schemes, and also a 2-D example using a 2"''-order SP-BGK scheme. Each of the examples is very sensitive to 
the accuracy of the scheme. Some of the tests run for millions of numerical steps. If the scheme is not a well-balanced 
one, the accumulation of any small numerical error would become significant for such a long time integration ifioll . 

6.1. Shock tube under gravitational field 

This case is the standard Sod test under gravitational field. The computational domain is x e [0, 1] which is divided 
into 100 cells. Reflection boundary condition is used on both ends. The initial condition is 

p= l.O,U ^0.0,p = 1.0forx<0.5, 

and 

p = 0.125, U = 0.0, p = 0.1 for x > 0.5. 

The gravitational force G takes a value G = - 1 .0 in the x-direction. So the potential jump at each cell interface 
becomes 

A0 = -GAx = 0.01. 

The computational results at f = 0.2 are presented in fig.|5] |6]for the density, pressure and velocity from the P'-order 
SP-KFVS, 1*' and 2"''-order SP-BGK schemes. From these figures, we can find that SP-KFVS scheme has larger 
numerical dissipation than that in SP-BGK scheme, and l*'-order scheme is more dissipative than 2"''-order one. The 
results calculated by the 2"'' order SP-BGK scheme fits the exact solution very well. Due to the gravitational force, the 
density distribution inside the tube is pulled back in the negative x-direction. In some region, the flow velocity even 
becomes negative. 
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6.2. Isolated gravitational system with adiabatic wall 

The second test case is also on a computational domain x e [0, 1] with 50 cells. There are limited number of 
gravitational potential jumps at locations x - 0.21, 0.41, 0.61 and 0.81 with a large value 

AO = 2.0. 

The initial flow distributions inside the domain has constant values of 

p = \.Q,pU = 0.0, and pE = 2.5. 

After a long time (f = 1000), the flow distributions settle down into a piecewise constant state which are shown in 
the fist picture of fig. |7J where the symbols are the numerical solutions and the solid lines are the exact hydrostatic 
solutions. The velocity distributions are also shown in fig. |2l For the 1" order schemes, the oscillation of velocity 
around zero is on the order of 10"^. This is mainly caused by the error in numerical integrations because there is no 
exact solution for most integrals in Ea.(l92]i-(l95Tl. In fact, the precision of numerical integration for the integrals is on 
the order of lO"*" -10"^. Since the potential jumps are large and the high order scheme uses more integral evaluations, 
the velocity distribution calculated by 2'"^ order scheme is a little bit worse than the Ist-order ones. If a better accuracy 
can be achieved for the numerical evaluation of the integrals, the velocity error can be further reduced to machine zero. 

6.3. Perturbation of the ID isothermal equilibrium solution 

This test case is from LeVeque and Bale's paper • We consider an ideal gas with y - 1 .4 on an initial isothermal 
hydrostatic state, 

poW = PoW = e''\ and Uo{x) = 0, 
for X e [0, 1]. Initially, the pressure is perturbed by 

p(x,t^Q) = poW + T/e^^^'"^"'', 

where a - 100, jco = 0.5 and rj is the amplitude of the perturbation. The gravitational field is the same as in example 
16.11 The computation is conducted with 100 grid points in the whole domain and stops at time t - 0.25. Fig. [8] show 
the results from SP-KFVS and SP-BGK schemes, where SP-KFVS has larger numerical dissipation than SP-BGK 
scheme. The results calculated by the 2'"^ order SP-BGK scheme matches the exact solution very well. 

Also in fig. |9] we show the convergency rate of our 2'"'-order SP-BGK scheme, where the number of cells is N and 
the error is the L°° error. From the figures, we can conclude our 2'"'-order SP-BGK scheme has a 2nd-order accuracy 
even with the modeling of piecewise constant potential. 

6.4. One-dimension gas falling into a fixed external potential. 

This case is taken from the paper by Slyz and Prendergast |@] to investigate the numerical accuracy of the BGK 
scheme. The gas is initially stationary (U - 0) and homogeneous {p - \, e - \, where e is the internal energy). The 
gravitational potential has the form of a sine wave, 

L . Inx 
= -00— sm— , 
2n L 

where L = 64 is the length of the computational domain and (pQ - 0.02. The ratio of the specific heat 7 = 5/3. The 
periodic boundary conditions are implemented in this system. Simulation results are presented with Ax = 1 and at the 
output time t - 250000 (more than 500000 time steps). After the initial transition, the system is expected to reach an 
isothermal hydrostatic distribution, where the temperature settles to a constant with zero velocity, i.e., 

T{x, t) = To, and f/ = 0. 

The velocity and temperature distributions computed by different symplecticity preserving schemes are shown in 
fig- E] [l2] The numerical error is smaller than that in 1 10]. Moreover, the results can be further improved if a better 
numerical integration for the integral evaluation can be adopted. 
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6.5. Rayleigh-Taylor instability. 

This test case also comes from 10] • Consider an isothermal equilibrium idea gas {y - 1 .4) in a 2D polar coordinate 

po(r) = poir) = l:^^-"*^^'''', f/o = 0, 

a 

where 

a = 2.68, ro = 0.258 forr < n, T n = 0.6(1 + 0.02 cos(206»)) foxdensity, 

and <{ 

a = 5.53, ro = -0.308 forr > n, [ n = 0.62324965 for pressure. 

The potential satisfies -V(p{r) - 1.5. The time evolutions of the density distributions at times t - 0,0.8, 1.4 and 2.0 
are shown in fig. [13] Fig. [T4]shows a scatter plot of the density as a function of the radius. These figures clearly show 
that the hydrostatic solution can be well kept and the flow motion is limited around the unstable interface. 

7. Conclusion 

In this paper, based on the Liouville's theorem and symplecticity-preserving property of a Hamiltonian flow, a 
wefl-balanced gas-kinetic BGK scheme (SP-BGK) has been developed for a hydrodynamic system under gravitational 
field with the modeling of piecewise constant potentials. As shown in the paper, in order to design such a scheme, the 
equilibrium state used has to be an exact Maxwellian distribution function. At the same time, the physical mechanism 
of particle transport across a potential barrier has to be explicitly followed in the equilibrium states modeling and 
the flux evaluation. As far as we know, the method presented in this paper is the first exact well-balanced scheme 
for the Navier-Stokes equations under gravitational field. At the same time, the particle transport mechanism across 
a potential jump in the current kinetic formulation follows the physical principles closely, which is valid under any 
general physical situation. Both the shock capturing and well-balanced properties are automatically obtained under the 
corresponding physical conditions. Mathematically, it has been proved that the SP-BGK method is a well-balanced 
scheme which could keep the hydrostatic state forever. In this paper, the design of the well-balanced scheme comes 
from the first principles of physics, instead of using the well-balanced condition as the starting point in the design of 
such a scheme. 
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Appendix 

Formulae in the two-dimensional case: 

1 . Equilibrium states 
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Case I. (/>j < (/>j+i, define Uc = ^j2{<pj+\ - (f)j). 



du dv d^ 



( 1 ^ 

-u 

V 



dwdvcJ^ 



+ ///-co /j+i 0' V, 



U 



u c2 



du dv d^. 



u 

m 



du dvd^ 



+ ///-I /7+l(-^7+l/2. 0. «. V, ^) 



dw dv d^. 



Case 2. > define Vc - y2((^y - 

^■.i/2 = //F//^;-i/2.0,«,^) 



l(„2 + ^2+^2) ; 



dwdvd^ 



+ ///-cr'/7+l(-^;+l/2,0,M,^) 



U 



u el 



du dv d^. 



(80) 



(81) 



(82) 
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^.l/2 = //F/X^7>l/2.0,«,^) 



u 

uv 



dwdvd^ 



/ 1 ^ 
-u 



dwdvd^ 



( 1 ^ 

u 



du dv d^. 



2. Fluxes 

Case I. 4>j < 4>j+i' define Uc = -^2{^j+\ - (pj). 

^kl/2« = //F/;(^.-.l/2,^,«,^) 



U \ 

uv 



dwdvd^ 



(83) 



—u 

-UV 



dwdvd^ 



(84) 



•///-lo/'+l(-^J+l/2'^'«'^) 



-U yju^ + 



Hu(u^ + U^) + uv^ + uf) ) 



^ 2 



dM dv d^. 



^;.i/2(o=//r;//^.-.i/2,^,«,^) 



UV 

{ Uuiu^ - U^) + uv^ + uf) ) 



dwdvd^ 



+ //il/7+l(-^J+l/2'f' 



U 
UV 

\{U^ +MV^ + M^^) ) 



du dv d^. 



(85) 
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Case 2. > </>j+i, define Uc = -y/2(0y - (f>j+i). 



u 

UV 



Au dv d^ 



+ ///-J'/7+l(^7+l/2,?,M,^) 



-U - f/2 

MV 



dM dv d^. 



UV 



dwdvd^ 



-M 

m2 

-MV 

k-M^ - MV^ - M^2) 



dwdvd^ 



+ //il/i+i(^J+i/2.^,«.^) 



MV 

j(m^ + MV^ + U^^) , 



dw dv d^. 



Formulae in the one-dimensional case: 

1. Equilibrium states: 

Case \. (pj< <pj+i, define Uc = ^2(0^+1 - (pj). 

W^j.l/2 = //o""//^.--l/2>0,«,^) 



^ 1 



+ //o'''//^7+l/2.0,M,^) 



V 2 



-M 



dwd^ 



V«2 + [/2 



dw d^. 
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u 



1 ^ 
U 



Case 2. 0^- > (pj+i, define Uc = ^j'2-{<pi - (f>j+i). 

^■.l/2 = /F//^.--l/2.0,M,^) 



( 1 



+ 1 L"' fMxM/2,0,u,O 



^lu^-Ul 



dw d^. 



dud^ 



r 



1 ^ 



dud^ 



+ 1 j^^fMxnU2,0,u,^) 



u 

w+e) ) 



du d^. 



2. Fluxes: 

Case 1. 0; < (pj+i, define Uc = ■\j2{<pj+\ - (f)j). 



dud^ 



+ JJo"'fjixM/2,t,u,0 



-u 



dwd^ 



V 2 



Hu{u^+u^)+ue) ) 



dwd^. 
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^ 2 



U - t/2 



dM d^. 



(93) 



Case 2. 0y > define Uc - ^%<Pj - <Pi+\)- 

^;.i/2«=/ir/.(^;.i/2,f,«,^) 



dM d^ 



+ //-»' /)+l(-*y+l/2'f'«'^) 



dMdf. 



(94) 



^;.l/2(0 = //o''"/.(^7.1/2,f,«,^) 



V 2 



M -y/M2 + [/2 
i(M(M2 + J/2) + Mf2) 



dwd^ 



+ //-V//+l^-^i+l/2'^'«'^) 



-M 

Ic .,3 .,e2 



^ 2 



(-«' - Mr) j 



dwd^ 



(95) 



+ /ll/>+l(-^^+l/2'f'«'f) 



i(M^ + <2) 



dwd^. 



Remarks on the integral evaluation: in the above formulae, there are many integrals which can not be analytically 
evaluated, e.g., , " )du. Therefore, a numerical integration method in 17] has been used. 
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Figure 1 : Reconstruction of the conservative variables at the cell interface. 




j y+1/2 7+1 

Figure 2: The modeling of the initial and equilibrium distribution functions at the cell interface for the BGK scheme without external forcing field. 



Figure 3: The particles' movement at the interface with potential jump 0^ < (f>j+\. 
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Figure 4: The modeling of the initial and equilibrium distribution functions at the cell interface for the SP-BGK scheme with a potential jump at 
the cell interface. 
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Figure 5: Density distributions for the shock tube problem under gravitational field. 
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Figure 6: Pressure and velocity distributions for the shock tube problem under gravitational field. 
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Figure 7: The first figure shows the Density, pressure and velocity distributions calculated by 2'"' order SP-BGK for isolated gravi- 
tational system with adiabatic wall. Other figures are velocity distributions in this test case. 
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Figure 8; Perturbation of pressure on an isothermal equilibrium solution. Left; = 0.01; right: rj = 0.001. 
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Figure 10: Density distribution calculated by 2"''-order SP-BGK for gas falling into a fixed external potential in 1-D case. 
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Figure 11: Velocity distributions for gas falling into a fixed external potential in 1-D case. The exact solution should have a zero 
velocity. The error is due to the numerical integration, e.g., I g(u)( , " )du, where there is no analytic solution. 
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Figure 12: Temperature distributions for gas falling into a fixed external potential in 1-D case. 



1 I- 




Figure 13: Rayleigh-Taylor instability with gravitational field directed radially inward. Density contours at time / = 0, 0.8, 1.4, 2.0 
are shown in the four quadrants, starting with the initial data in the upper right corner and progressing clockwise. 
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